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ABSTRACT 

A series system is simulated to obtain lower confidence 
limits on system reliability using Bayesian techniques. A 
comparison between classical and Bayesian methods is made. 
Random beta variate generators are developed and used in the 
simulation. The results of the simulation are tabulated for 
easy comparison of the Bayesian and classical methods. The 
values of lower confidence limits that are realized using 
the Bayesian method decrease as the number of components 
increase. In most cases, as the number of components 
increase, the Bayesian method appears to yield lower values 
of loiter confidence limits than the classical method. 
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I. INTRODUCTION 



In a series system the reliability of the system, RS, 
depends upon the reliability of each of the k components or 
subsystems of which it is composed. This can be modeled as 
follows : 

k 

RS = ir R. , 
i = l 1 

where R^ is the reliability of the ith component. Using 
success/fail test data classical or Bayesian methods can be 
used to compute lower confidence limits on RS. In a simple 
system with few components the Bayesian method offers many 
advantages over the more conservative classical method. 
However, as the system becomes more complex with a higher 
number of components, problems start to arise with the use 
of standard Bayesian techniques. 

The purpose of this paper is to compare classical and 
Bayesian methods of computing lower confidence limits. Com- 
puter simulations were used to determine lower confidence 
limits on system reliability using Bayesian methods. A 
Poisson approximation was used to obtain lower confidence 
limits for the classical approach. Several factors were 
varied but the main interest of the investigation was the 
effect of the number of components of a system on the values 
of lower confidence limits on system reliability. 
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II. METHODS 



In order to compute lower confidence limits on system 
reliability of a series system two methods can be used. 

The first method to be discussed is the classical approach 
which is based solely on test results. The second method 
is the Bayesian approach which utilizes assumptions based 
on prior knowledge of similar systems. 

A. CLASSICAL METHOD 

The classical approach to system reliability uses only 
the mission test data on each component. These results are 
then combined in order to obtain an estimate of the system 
reliability. In the most simple case where k components are 
each tested n times and there are no failures in any of the 
components it is assumed that this is equivalent to testing 
the entire system n times with no failures. In this case 
the lower 100(l-a)% confidence limit can be found in the 
same manner as with just one item. This is done by solving 
for p in the equation: 

E (hphl-p) 1 ' ■ a 

j-s 1 

where s = number of successes. 

The solution is a lower 100(l-a)l confidence limit on RS. 

N 

In the case of zero failures RS T , = /a . When only one 

L (a) 

failure occurs the same procedure can be used. However, 
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when more than one failure occurs the procedure becomes 
more complicated. 

In the case where each component is tested n times and 
there are few component failures, an approximation to the 
classical value of RS T , s can be found as follows. 



Let = 1-R^ , 



k 

then RS = Z (1-Q. ) . 

i=l 1 



If the Q^’s are small then 



k 

RS = 1 - Z Q. . 

• -i 1 



Let f^ = the number of failures on the ith component. 
Since f. is a sum of n Bernoulli trials it is binomial 

l 

(n,Q^). If CK is small then each f^ is approximately 
Poisson (nO^) . 



k 

Let F = Z f- ; 
i-1 1 

k 

then F is approximately Poisson (n Z Q.)> 

i = l 1 

and the upper lOO(l-a)! confidence limit [Ref. 1] for 
k X*, 2 (F+l) 

Z Q. is = , where y 2 » 2(F + L) is the 100(l-a)th 

^ i / n ot 

percentile of the x 2 distribution with 2 (F + l) degrees of 
freedom. 
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X*, 2 (F+l) 

Thus RS L(a) - 1 - 

is an approximation to the classical value for a lower 
100(l-cc)% confidence limit on RS. 

B. BAYESIAN METHOD 

The Bayesian approach to system reliability is to treat 
each component reliability, R^, as a random variable. One 
particular method is to assume that each component relia- 
bility has a beta distribution. The beta density function 
can be defined as follows: 

F(r;a,b) = BTaTFT r a " 1 (l-r) b " 1 for 0 < r < 1 , a >0 , b > 0, 

where B(a,b) = • 

A prior beta density B(r;a,b) is assumed for each component 
before testing begins. The test data is then used to obtain 
a posterior beta density B(r;a+s,b+f) where s is the number 
of successful tests and f is the number of failures for the 
component . 

By generating a random beta variate for each component 
based on its posterior 'density a value of RS is calculated 
using the model : 

k 

RS = Z R. . 
i = l 1 

A distribution of RS values is simulated by repeating this 
procedure. The lower confidence limits are thus estimated 
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by the appropriate percentile points of the simulated 
distributions of RS values. 

The main advantage of the Bayesian method is that past 
experience with technically similar hardware can be used to 
determine appropriate priors for each component. This allows 
higher reliability goals to be met with less testing. How- 
ever, some precautions must be observed in selecting the 
beta priors. Choosing a beta prior of B(r;a,b), where a and 
b are integers, is equivalent in the classical sense to 
assuming that (a+b-1) tests have been performed and that a 
of these were successes. It is equivalent because the ath 
percentile point of the B(r;a,b) distribution is the lower 
lOO(l-ct)! confidence limit for p when there are (a+b-1) 
Bernoulli trials, each with probability p of success, and a 
of these trials are successful. This is readily apparent 
from an identity that relates the beta distribution to the 
binomial distribution [Ref. 2]. 

As the sum of the beta parameters (a+b) in the prior 
density increases relative to the number of tests to be 
observed the resulting lower confidence limits are more of 
a function of the assumed prior density than of the test 
results. One procedure is to limit the sum of the two 
parameters to a function of the number of tests. If the b 
parameter is less than one a different problem exists even 
though the sum of the parameters is small. This problem 
can be shown by the following discussion. 
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For any a, a family of beta distributions exist such 
that the ath percentile point of each is the same. If the 
b parameter is set equal to 1, then for any B(r;a,b) 
distribution an a* can be found such that the ath percen- 
tile point of B(r;a,b) and B(r;a*,l) are the same. The 
B(r;a*,l) distribution is convenient since it can be related 
to the classical case of having already tested a* items with 
no failures. The problem when the b parameter is less than 
one can be shown by converting the initial B(r;a,b) prior to 
a B(r;a*,l) prior. In this case if there are no failures 
the posterior distribution will be B (r ; a* + s , 1) where s is 
the number of successful tests. If the initial beta prior 
is not converted the posterior distribution will be 
B(r;a+s,b), which can be converted to B (r ; (a+s) * , 1) . In 
all cases that were investigated the (a+s)* was larger than 
a*+s when the b parameter was less than one. Thus, in an 
equivalent classical sense, it appears that each successful 
test will be counted as more than one success. This 
amplification effect seems to increase rapidly as the b 
parameter approaches zero. Also if the beta prior (B(r;a,b), 
where b is less than one, is converted to a beta B(r;a*,l) 
to ensure that each success counts only as one success, the 
resulting a* may be too large. The details of the conversion 
are given in the next section and a table of the resulting 
a*'s for various values of b can be found in Appendix A. 
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III. SIMULATION 



The purpose of this investigation was to compare lower 
confidence limits derived by using Bayesian methods with 
those obtained by using the classical method. In order to 
simplify the simulation and to be able to make comparisons 
it was assumed that each of the components had the same 
beta prior density and that each component was mission 
tested the same number of times. Using the same beta prior 
for each component is not a necessary condition for the 
simulation to work. 

A. SIMULATION PROCEDURE 

A series system of k components was simulated in the 
following manner. Let Bl^(r;a^,b^) be the initial beta 
prior for R^ . For this investigation the following three 
beta priors were used: B (r ; 5 . 0 , 0 . 03) , B (r ; 5 . 0 , 0 . 05) and 

B (r ; 5 . 0 , 0 . 1 0) . For each Bl^(r;a^,b^) a new beta prior, 
B2^(r;a^,l), was computed such that the 10th percentile 
point of B2^(r;a^,l) was the same as the 10th percentile 
point of B1 ^ (r ; a^ , b ^) . The value for a^ is easily computed 
since if the b parameter is 1 then the percentile point (P) 
of a B(x;a,l) distribution is given by the equation P = x . 
Thus a? = ln(0.10)/ln x, where x is the 10th percentile 
point of B1 . (r ; a - , b . ) . 

Let n^ = the number of trials for each component, 
and s^ = the number of successes for each component. 
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The B2^(r;a*,l) prior was then adjusted in the following 
manner : 

B21^(r;a^,l) where a^ = min (0 . 75n^ , a^) 

B22^(r;a^,l) where a^ = min(1.0 n^,a*) 

B23^(r;a^,l) where a^ = min(1.5 n^,a^) . 

A beta posterior B3^(r;a^,b^) for each was computed in 
four cases as follows: 



Case 


1 


B3 i (r;a i 


•V 


= B1 . (r ; a . +s • , b . +n . - s . ) 

i v ’i i’i l i J 


Case 


2 


B3 i (r;a i 


■V 


= B21 . (r;a . +s . , 1+n . - s . ) 
1 v ’ 1 1 ’ 1 x J 


Case 


3 


B3 i (r;a . 


>v 


= B22.(r; a.+s.,l+n.-s.) 
r ’ i i ’ i i J 


Case 


4 


B3 i (r;a i 


.bp 


= B23 . (r ; a . +s . , 1+n . - s . ) 

i v > i 1 > 1 



By generating B3^(r;a^,b^) random variates, which represent 
the posterior distributions of each of the components, a 
value for RS was obtained by taking their product. This 
procedure was repeated 500 times and realized values for RS 
were then ordered. The lower 100(l-a)% confidence limits 
were then determined by selecting the appropriate percentile 
points of the simulated distribution on RS. 

B. GENERATION OF RANDOM BETA VARIATES 

Two means of generating random beta variates were used 
in the simulation. In the case where the b parameter of the 
B(r;a,b) distribution was an integer an exponential genera- 
tor was used to realize random beta variates. In the case 
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of noninteger parameters, Monte Carlo rejection techniques 
were used. 

In order to generate random B(r;a,b) variates where b 
is an integer the following logic was used [Ref. 3] . Assume 
Y is B(y;a,b) . 

Define U = -InY. 

Then the moment generating function 

My t = E [e tU ] = E[e" lnY ] = E[Y _t ] . 

Thus 

p r y * i i 1 a — t — 1 / -i \ b “ 1 i 

E < Y ] = BTT7BT / y (1 ' y) dy > 

B (a- 1 ,b) _ r (a+b) . r(b)r(a-t) 

" B (a , b") rCaj'-T'CbT T(a + b-tJ 5 

a (a + 1 ) (a+b- 1) 

(a-tj (a+1- t} . . . . (a+b-l-tj 

a+b-1 , 



Since 

I 

is the moment generating function for an exponential random 
variable with parameter X, then U is the sum of b indepen- 
dent exponential random variables with failure rates a, 
a+1, ..., a+b-1, 

b-1 

i.e., U = £ U. where U. is exp (a+j). 

j=0 - 1 J 
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Since Y = e random beta variates can be generated with 
the use of existing exponential generators. For the simu- 
lation, where the posterior distributions on were 
B3 • fr ; a . ,b . ) and the b.'s were integers, values for RS 
were obtained by using the following formula: 

k V 1 

-InRS = £ E U. . , 

i = l j = 0 1J 

where 

IKj is exponential (a^ + j). 

In the case where the b parameter of the beta density 
was noninteger, random beta variates were generated by 
using Monte Carlo rejection techniques [Ref,. 4], Since 
this investigation was concerned with small values of the b 
parameter the Monte Carlo technique had to be modified 
slightly. The procedure that was used is outlined below. 



1. Compute P 

2. Set H 

3. Generate R 

4. If R > P 

5. Generate R1 

6. If R1 > 0.999 

7. Generate R2 

8. If (R2-H < f x (Rl) 

9. Go to 5. 



P=Prob(X £ 0.999) where X is B(x;a,b). 
H=Max value of f x (x) f° r 0 £ x £ 0.999. 
R=Random U(0,1) number. 

Then X-=1.0, Go to 3. 

Rl=Random U(0,1) number. 

Then Go to 5. 

R2=Random U(0,1) number. 

Then X.=R1, Go to 3. 

l ’ 



This procedure was repeated until 500 random beta variates 
were realized for each posterior density. 
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C. MODIFICATION OF SUBROUTINE BDTR 



Subroutine BDTR [Ref. 5] was used in the simulation to 
convert beta priors and to compute beta densities in the 
generation of random beta variates using Monte Carlo tech- 
niques. The subroutine BDTR computes the probability that 
the random variable U is less than or equal to x, where U 
is distributed according to the beta (u;a,b) distribution 
and 0 _< x £ 1 . The value of the density at x is also com- 
puted. In order for the computations to be valid the sum 
of the parameters must be greater than or equal to one. 

In the subroutine BDTR this condition is met by restricting 
both parameters to values greater than or equal to 0.5. 

The subroutine was modified so that smaller values of the 
b parameter could be used since the a parameter was always 
greater than or equal to one. Two additional variables 
were added to the parameter list of the subroutine which 
increased the efficiency of the simulation by acting as 
flags’ and saving values when the subroutine was being 
called during the generation of random beta variates. All 
modifications to the subroutine are shown in the listing 
of the computer program. 
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IV. RESULTS 



The results of the simulation are given in tables for 
each of the three beta priors that were used. The number 
of failures was the total number of failures in all of the 
components. For n=50, and in a similar manner where n=75 
and n=100, the failures were assigned as follows: 



Failures 


Successes 






0 


s i 


= 


50, 


i 1,2,. 


..,k 


1 


s . 
1 


= 


49, 

50, 


i=l 

i=2 , 3 , . 


. . ,k 


2 


S, 


= 


49, 


i=l , 2 






1 


= 


50, 


i = 3, 4, . 


. . , k 


3 


s . 


= 


48, 


i = l 






1 


= 


49, 


i = 2 








= 


50, 


i — 3 , 4 , . 


..,k 


4 


s . 
1 


= 


48, 

50, 


i = l , 2 
i=3,4, . 


. . , k 


5 


s . 
1 


= 


48, 

49, 


i = l , 2 
i = 3 








= 


50, 


• 

#V 

LO 

II 

•H 


..,k 


6 


S • 


= 


48, 


i = l , 2 






1 


= 


49, 


i = 3 , 4 








= 


50, 


i= 5 , 6 , . 


..,k 



The beta priors for the Bayesian cases were as follows: 



Case 1 Bl(r;a,b) 

Case 2 B21(r;a,l) a 

Case 3 B22(r;a,l) a 

B23(r;a,l) a 



min (0 . 7 5n , a* ) 
min (1.0 n,a*) 
min(l . 5 n,a*) 



n=50, 75,100 
n=50, 75,100 
n=50, 75,100 



Case 4 
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The appropriate values for Bl(r;a,b) and B2(r;a*,l) are 
given at the top of each table. 
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90% Lower Confidence Limits on RS 



Bl(r; 


5.0,0. 


03) converts 


to B2(r; 


596.8,1) 






Fifty 


tests 


for each component 












Classical 




Bayes 


ian 




k Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


40 


0 


0.9540 


0.9532 


0.5751 


0.6162 


0.6789 


40 


1 


0.9222 


0.9272 


0.5739 


0.6151 


0.6779 


40 


2 


0.8936 


0.9038 


0.5599 


0.6020 


0.6663 


40 


3 


0.8664 


0.8827 


0.5541 


0.5966 


0.6116 


40 


4 


0.8401 


0.8647 


0.5493 


0.5920 


0.6575 


40 


5 


0.8145 


0.8410 


0.5429 


0.5861 


0.6523 


40 


6 


0.7894 


0.8247 


0.5338 


0.5775 


0.6446 


30 


0 


0.9540 


0.9616 


0.6547 


0.6903 


0.7434 


30 


1 


0.9222 


0.9345 


0.6467 


0.6829 


0.7371 


30 


2 


0.8936 


0.9119 


0.6384 


0.6752 


0.7305 


30 


3 


0.8664 


0.8897 


0.6258 


0.6636 


0.7205 


30 


4 


0.8401 


0.8696 


0.6247 


0.6627 


0.7196 


30 


5 


0.8145 


0.8491 


0.6179 


0.6563 


0.7140 


30 


6 


0.7894 


0.8307 


0.6075 


0.6467 


0.7057 


20 


0 


0.9540 


0.9731 


0.7452 


0.7731 


0.8139 


20 


1 


0.9222 


0.9398 


0.7366 


0.7653 


0.8074 


20 


2 


0.8936 


0.9159 


0.7277 


0.7572 


0.8006 


20 


3 


0.8664 


0.8950 


0.7192 


0.7495 


0.7940 


20 


4 


0.8401 


0.8754 


0.7074 


0.7388 


0.7850 


20 


5 


0.8145 


0.8549 


0.7012 


0.7331 


0.7802 


20 


6 


0.7894 


0.8348 


0.6924 


0.7251 


0.7734 
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90% Lower Confidence Limits on RS 



B1 (r ; 5 


.0,0.03) 


converts 


to B2(r; 


596.8,1) 






Fifty 


tests for 


each component 










Classical 




Bayes ian 




k Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


10 


0 


0.9540 


0.9840 


0.8478 


0.8654 


0.8908 


10 


1 


0.9222 


0.9474 


0.8413 


0.8597 


0.8861 


10 


2 


0.8936 


0.9217 


0.8306 


0.8501 


0.8782 


10 


3 


0.8664 


0.9023 


0.8156 


0.8368 


0.8673 


10 


4 


0.8401 


0.8817 


0.8069 


0.8289 


0.8607 


10 


5 


0.8145 


0.8576 


0.7949 


0.8182 


0.8538 


10 


6 


0.7894 


0.8399 


0.7822 


0.8067 


0.8423 


5 


0 


0.9540 


0.9924 


0.9163 


0.9264 


0.9407 


5 


1 


0.9222 


0.9521 


0.9023 


0.9140 


0.9306 


5 


2 


0.8936 


0.9259 


0.8850 


0.8986 


0.9180 


5 


3 


0.8664 


0.9062 


0.8763 


0.8910 


0.9119 


5 


4 


0.8401 


0.8833 


0.8637 


0.8798 


0.9027 


5 


5 


0.8145 


0.8616 


0.8466 


0.8645 


0.8902 


5 


6 


0.7894 


0.8435 


0.8400 


0.8587 


0.8856 
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90% Lower Confidence Limits on RS 



B1 (r ; 5 . 0 , 0 . 03) converts to B2 (r ; 596 . 8 , 1) 
Seventy-five tests for each component 

Classical Bayesian 



k 


Failures 


Case 


Case 1 


Case 2 


Case 3 


Cas e 4 


40 


0 


0.9693 


0.9701 


0.6915 


0.7242 


0.7724 


40 


1 


0.9481 


0.9496 


0.6906 


0.7233 


0.7717 


40 


2 


0.9290 


0.9366 


0.6794 


0.7130 


0.7629 


40 


3 


0.9109 


0.9135 


0.6747 


0.7088 


0.7593 


40 


4 


0.8934 


0.9032 


0.6708 


0.7051 


0.7561 


40 


5 


0.8763 


0.8865 


0.6658 


0.7006 


0.7523 


40 


6 


0.8596 


0.8721 


0.6582 


0.6936 


0.7463 


30 


0 


0.9693 


0.9756 


0.7540 


0.7811 


0.8206 


30 


1 


0.9481 


0.9552 


0.7478 


0.7755 


0.8160 


30 


2 


0.9290 


0.9402 


0.7415 


0.7697 


0.8111 


30 


3 


0.9109 


0.9186 


0.7318 


0.7610 


0.8038 


30 


4 


0.8934 


0.9062 


0.7310 


0.7602 


0.8031 


30 


5 


0.8763 


0.8920 


0.7256 


0.7553 


0.7989 


30 


6 


0.8596 


0.8780 


0.7175 


0.7480 


0.7927 


20 


0 


0.9693 


0.9803 


0.8219 


0.8423 


0.8717 


20 


1 


0.9481 


0.9591 


0.8157 


0.8367 


0.8671 


20 


2 


0.9290 


0.9435 


0.8091 


0.8309 


0.8622 


20 


3 


0.9109 


0.9235 


0.8028 


0.8251 


0.8575 


20 


4 


0.8934 


0.9093 


0.7941 


0.8174 


0.8510 


20 


5 


0.8763 


0.8958 


0.7895 


0.8132 


0.8476 


20 


6 


0.8596 


0.8804 


0.7829 


0.8073 


0.8427 
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90% Lower Confidence Limits on RS 



B1 (r ; 5 . 0 , 0 . 03) converts to B2 (r ; 596 . 8 , 1) 
Seventy- five tests for each component 

Classical Bayesian 



k 


Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


10 


0 


0.9693 


0.9891 


0.8957 


0.9082 


0.9258 


10 


1 


0.9481 


0.9647 


0.8912 


0.9041 


0.9225 


10 


2 


0.9290 


0.9474 


0.8836 


0.8974 


0.9171 


10 


3 


0.9109 


0.9267 


0.8732 


0.8882 


0.9095 


10 


4 


0.8934 


0.9138 


0.8669 


0.8825 


0.9049 


10 


5 


0.8763 


0.9001 


0.8584 


0.8750 


0.8987 


10 


6 


0.8596 


0.8836 


0.8492 


0.8668 


0.8920 


5 


0 


0.9693 


0.9940 


0.9434 


0.9503 


0.9600 


5 


1 


0.9481 


0.9677 


0.9338 


0.9418 


0.9532 


5 


2 


0.9290 


0.9504 


0.9218 


0.9313 


0.9447 


5 


3 


0.9109 


0.9304 


0.9160 


0.9261 


0.9405 


5 


4 


0.8934 


0.9159 


0.9072 


0.9183 


0.9342 


5 


5 


0.8763 


0.9002 


0.8951 


0.9077 


0.9255 


5 


6 


0.8596 


0.8850 


0.8908 


0.9038 


0.9224 
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901 Lower Confidence Limits on RS 



B1 (r ; 5 . 0 , 0 . 03) converts to B2 (r ; 596 . 8 , 1) 

One hundred tests for each component 

Classical Bayesian 



k 


Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


40 


0 


0.9770 


0.9749 


0.7583 


0.7850 


0.8239 


40 


1 


0 . 9611 


0.9586 


0.7576 


0.7843 


0.8234 


40 


2 


0.9458 


0.9481 


0.7483 


0.7759 


0.8163 


40 


3 


0.9332 


0.9380 


0.7445 


0.7725 


0.8134 


40 


4 


0.9200 


0.9266 


0.7412 


0.7695 


0.8109 


40 


5 


0.9073 


0.9106 


0.7372 


0.7658 


0.8079 


40 


6 


0.8947 


0.9019 


0.7308 


0.7601 


0.8030 


30 


0 


0.9770 


0.9817 


0.8091 


0.8308 


0.8622 


30 


1 


0.9611 


0.9641 


0.8042 


0.8264 


0.8585 


30 


2 


0.9458 


0.9514 


0.7991 


0.8218 


0.8547 


30 


3 


0.9332 


0.9417 


0.7913 


0.8149 


0.8490 


30 


4 


0.9200 


0.9270 


0.7906 


0.8142 


0.8484 


30 


5 


0.9073 


0.9133 


0.7862 


0.8102 


0.8451 


30 


6 


0.8947 


0.9070 


0.7797 


0.8043 


0.8402 


20 


0 


0.9770 


0.9861 


0.8632 


0.8793 


0.9022 


20 


1 


0.9611 


0.9685 


0.8583 


0.8749 


0.8986 


20 


2 


0.9458 


0.9541 


0.8532 


0.8703 


0.8948 


20 


3 


0.9332 


0.9446 


0.8481 


0.9658 


0.8911 


20 


4 


0.9200 


0.9304 


0.8413 


0.8596 


0.8861 


20 


5 


0.9073 


0.9183 


0.8337 


0.8564 


0.8834 


20 


6 


0.8947 


0.9090 


0.8324 


0.8518 


0.8796 
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901 Lower Confidence Limits on RS 



B1 (r ; 5 . 0 , 0 . 03) converts to B2 (r ; 596 . 8 , 1) 
One hundred tests for each component 

Classical Bayesian 



k 


Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


10 


0 


0.9770 


0.9918 


0.9207 


0.9303 


0.9438 


10 


1 


0.9611 


0.9727 


0.9172 


0.9272 


0.9413 


10 


2 


0.9458 


0.9595 


0.9114 


0.9220 


0.9371 


10 


3 


0.9332 


0.9474 


0.9034 


0.9150 


0.9314 


10 


4 


0.9200 


0.9336 


0.8984 


0.9106 


0.9278 


10 


5 


0.9073 


0.9210 


0.8919 


0.9047 


0.9127 


10 


6 


0.8947 


0.9101 


0.8847 


0.8984 


0.9179 


5 


0 


0.9770 


0.9956 


0.9572 


0.9625 


0.9699 


5 


1 


0.9611 


0.9753 


0.9499 


0.9560 


0.9647 


5 


2 


0.9458 


0.9609 


0.9408 


0.9480 


0.9582 


5 


3 


0.9332 


0.9504 


0.9364 


0.9441 


0.9551 


5 


4 


0.9200 


0.9344 


0.9296 


0.9382 


0.9503 


5 


5 


0.9073 


0.9222 


0.9204 


0.9300 


0.9436 


5 


6 


0.8947 


0.9121 


0.9171 


0.9271 


0.9413 
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901 Lower Confidence Limits on RS 



Bl(r; 


5.0,0.05) 


converts 


to B2 (r ; 


136.8,1) 






Fifty 


tests for 


each component 










Classical 




Bayesian 




k 


Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


40 


0 


0.9540 


0.9344 


0.5751 


0.6162 


0.6789 


40 


1 


0.9222 


0.9085 


0.5739 


0.6151 


0.6779 


40 


2 


0.8936 


0.8844 


0.5599 


0.6020 


0.6663 


40 


3 


0.8664 


0.8715 


0.5541 


0.5966 


0.6116 


40 


4 


0.8401 


0.8459 


0.5493 


0.5920 


0.6575 


40 


5 


0.8145 


0.8265 


0.5429 


0.5861 


0.6523 


40 


6 


0.7894 


0.8091 


0.5338 


0.5775 


0.6446 


30 


0 


0.9540 


0.9441 


0.6547 


0.6903 


0.7434 


30 


1 


0.9222 


0.9203 


0.6467 


0.6829 


0.7371 


30 


2 


0.8936 


0.8918 


0.6384 


0.6752 


0.7305 


30 


3 


0.8664 


0.8758 


0.6258 


0.6636 


0.7205 


30 


4 


0.8401 


0.8542 


0.6247 


0.6627 


0.7196 


30 


5 


0.8145 


0.8318 


0.6179 


0.6563 


0.7140 


30 


6 


0.7894 


0.8153 


0.6075 


0.6467 


0.7057 


20 


0 


0.9540 


0.9608 


0.7452 


0.7731 


0.8139 


20 


1 


0.9222 


0.9329 


0.7366 


0.7653 


0.8074 


20 


2 


0.8936 


0.9072 


0.7277 


0.7572 


0.8006 


20 


3 


0.8664 


0.8874 


0.7192 


0.7495 


0.7940 


20 


4 


0.8401 


0.8607 


0.7074 


0.7388 


0.7850 


20 


5 


0.8145 


0 . 8385 


0.7012 


0.7331 


0.7802 


20 


6 


0.7894 


0.8243 


0.6924 


0.7251 


0.7734 
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90% Lower Confidence Limits on RS 



Bl(r; 


5.0,0.05) 


converts 


to B2(r; 


136.8,1) 






Fifty 


tests for 


each component 










Classical 




Bayesian 




k 


Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


10 


0 


0.9540 


0.9749 


0.8478 


0.8654 


0.8908 


10 


1 


0.9222 


0.9452 


0.8413 


0.8597 


0.8861 


10 


2 


0.8936 


0.9142 


0.8306 


0.8501 


0.8782 


10 


3 


0.8664 


0.8934 


0.8156 


0.8368 


0.8673 


10 


4 


0.8401 


0.8672 


0.8069 


0.8289 


0.8607 


10 


5 


0.8145 


0.8469 


0.7949 


0.8182 


0.8538 


10 


6 


0.7894 


0 .8313 


0.7822 


0.8067 


0.8423 


5 


0 


0.9540 


0.9866 


0.9163 


0.9264 


0.9407 


5 


1 


0.9222 


0.9521 


0.9023 


0.9140 


0.9306 


5 


2 


0.8936 


0.9234 


0.8850 


0.8986 


0.9180 


5 


3 


0.8664 


0.9027 


0.8763 


0.8910 


0.9119 


5 


4 


0.8401 


0.8725 


0.8637 


0.8798 


0.9027 


5 


6 


0.7894 


0.8384 


0.8400 


0.8587 


0.8856 
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901 Lower Confidence Limits on RS 



B1 (r ; 5 . 0 , 0 . 05) converts to B2 (r ; 136 . 8 , 1) 
Seventy-five tests for each component 

Classical Bayesian 



k 


Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


40 


0 


0.9693 


0 . 9533 


0.6915 


0.7242 


0.7724 


40 


1 


0.9481 


0.9388 


0.6906 


0.7233 


0.7717 


40 


2 


0.9290 


0.9209 


0.6794 


0.7130 


0.7629 


40 


3 


0.9109 


0.9046 


0.6747 


0.7088 


0.7593 


40 


4 


0.8934 


0.8943 


0.6708 


0.7051 


0.7561 


40 


5 


0.8763 


0.8769 


0.6658 


0.7006 


0.7523 


40 


6 


0.8596 


0.8601 


0.6582 


0.6936 


0.7463 


30 


0 


0.9693 


0.9615 


0.7540 


0.7811 


0.8206 


30 


1 


0.9481 


0.9449 


0.7478 


0.7755 


0.8160 


30 


2 


0.9290 


0.9310 


0.7415 


0.7697 


0.8111 


30 


3 


0.9109 


0.9099 


0.7318 


0.7610 


0.8038 


30 


4 


0 . 8934 


0.8988 


0.7310 


0.7602 


0.8031 


30 


5 


0.8763 


0.8808 


0.7256 


0.7553 


0.7989 


30 


6 


0.8596 


0.8684 


0.7175 


0.7480 


0.7927 


20 


0 


0.9693 


0.9710 


„ 0.8219 


0.8423 


0.8717 


20 


1 


0.9481 


0.9526 


0.8157 


0.8367 


0.8671 


20 


2 


0.9290 


0.9367 


0.8091 


0.8309 


0.8622 


20 


3 


0.9109 


0.9172 


0.8028 


0.8251 


0.8575 


20 


4 


0.8934 


0.9057 


0.7941 


0.8174 


0.8510 


20 


5 


0.8763 


0.8885 


0.7895 


0.8132 


0.8476 


20 


6 


0.8596 


0.8770 


0.7829 


0.8073 


0.8427 
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90% Lower Confidence Limits on RS 



B1 (r ; 5 . 0 , 0 . 05) converts to B2 (r ; 136 . 8 , 1) 
Seventy-five tests for each component 

Classical Bayesian 



k 


Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


10 


0 


0.9693 


0.9826 


0.8957 


0.9082 


0.9258 


10 


1 


0.9481 


0.9597 


0.8912 


0.9041 


0.9225 


10 


2 


0.9290 


0.9421 


0.8836 


0.8974 


0.9171 


10 


3 


0.9109 


0.9245 


0.8732 


0.8882 


0.9095 


10 


4 


0.8934 


0.9101 


0.8669 


0.8825 


0.9049 


10 


5 


0.8763 


0.8959 


0.8584 


0.8750 


0.8987 


10 


6 


0.8596 


0.8819 


0.8492 


0.8668 


0.8920 


5 


0 


0.9693 


0.9892 


0.9431 


0.9503 


0.9600 


5 


1 


0.9481 


0.9657 


0.9338 


0.9418 


0.9532 


5 


2 


0.9290 


0.9476 


0.9218 


0.9313 


0.9447 


5 


3 


0.9109 


0.9286 


0.9160 


0.9261 


0.9405 


5 


4 


0.8934 


0.9160 


0.9072 


0.9183 


0.9342 


5 


5 


0.8763 


0.8972 


0.8951 


0.9077 


0.9255 


5 


6 


0.8596 


0.8842 


0.8908 


0.9038 


0.9224 
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901 Lower Confidence Limits on RS 



B1 (r ; 5 . 0 , 0 . 05) converts to B2 (r ; 1 36 . 8 , 1) 
One hundred tests for each component 

Classical Bayesian 



k 


Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


40 


0 


0.9770 


0.9649 


0.7583 


0.7850 


0.8151 


40 


1 


0.9611 


0.9507 


0.7576 


0.7843 


0.8145 


40 


2 


0.9458 


0 . 9373 


0.7483 


0.7759 


0.8071 


40 


3 


0.9332 


0.9277 


0.7445 


0.7695 


0.8015 


40 


4 


0.9200 


0.9148 


0.7412 


0.7695 


0.8015 


40 


5 


0.9073 


0.9086 


0.7372 


0.7658 


0.7983 


40 


6 


0.8947 


0.8945 


0.7308 


0.7691 


0.7932 


30 


0 


0.9770 


0.9721 


0.8091 


0.8391 


0.8308 


30 


1 


0.9611 


0.9579 


0.8042 


0.8264 


0.8513 


30 


2 


0.9458 


0.9444 


0.7991 


0.8218 


0.8473 


30 


3 


0.9332 


0.9336 


0.7913 


0.8149 


0.8413 


30 


4 


0.9200 


0.9207 


0.7906 


0.8142 


0.8406 


30 


5 


0.9073 


0.9120 


0.7862 


0.8102 


0.8372 


30 


6 


0.8947 


0.9011 


0.7797 


0.8043 


0.8313 


20 


0 


0.9770 


0.9794 


0.8632 


0.8793 


0.8970 


20 


1 


0.9611 


0.9641 


0.8583 


0.8749 


0.8932 


20 


2 


0.9458 


0.9496 


0.8532 


0.8703 


0 . 8893 


20 


3 


0.9332 


0.9395 


0.8481 


0.8658 


0.8854 


20 


4 


0.9200 


0.9270 


0.8413 


0.8596 


0 . 8801 


20 


5 


0.9073 


0.9172 


0.8377 


0.8564 


0.8774 


20 


6 


0.8947 


0.9063 


0.8324 


0.8518 


0.8733 
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90% Lower Confidence Limits on RS 



B1 (r ; 5 . 0 , 0 . 05) converts to B2 (r ; 136 . 8 , 1) 
One hundred tests for each component 

Classical Bayesian 



k 


Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


10 


0 


0.9770 


0.9871 


0.9207 


0.9303 


0.9408 


10 


1 


0.9611 


0.9696 


0.9172 


0.9272 


0.9382 


10 


2 


0.9458 


0.9566 


0.9114 


0.9220 


0.9338 


10 


3 


0.9332 


0.9457 


0 . 9034 


0.9150 


0.9277 


10 


4 


0.9200 


0.9327 


0.8984 


0.9106 


0.9240 


10 


5 


0.9073 


0.9219 


0.8919 


0.9047 


0.9190 


10 


6 


0.8947 


0.9091 


0.8847 


0.8984 


0.9135 


5 


0 


0.9770 


0.9935 


0.9572 


0.9625 


0.9682 


5 


1 


0.9611 


0.9739 


0.9499 


0.9560 


0.9627 


5 


2 


0.9458 


0.9599 


0.9408 


0.9480 


0.9560 


5 


3 


0.9332 


0.9478 


0.9364 


0.9441 


0.9526 


5 


4 


0.9200 


0.9359 


0.9296 


0.9382 


0.9476 


5 


5 


0.9073 


0.9271 


0.9204 


0.9300 


0.9406 


5 


6 


0.8947 


0.9114 


0.9171 


0.9271 


0.9381 
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90% Lower Confidence Limits on RS 



Bl(r; 


5.0,0.10) 


converts 


to B2 (r 


; 39 . 5 , 1) 






Fifty 


tests for 


each component 










Class ical 




Bayesian 




k 


Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


40 


0 


0.9540 


0.8825 


0.5751 


0.5821 


Same 














as 


40 


1 


0.9222 


0.8639 


0.5739 


0.5809 


Case 3 


40 


2 


0.8936 


0.8407 


0.5599 


0.5671 




40 


3 


0.8664 


0.8230 


0.5541 


0.5613 




40 


4 


0.8401 


0.8028 


0.5493 


0.5566 




40 


5 


0.8145 


0.7872 


0.5429 


0.5503 




40 


6 


0.7894 


0.7666 


0.5338 


0.5412 




30 


0 


0.9540 


0.9055 


0.6547 


0.6608 




30 


1 


0.9222 


0.8842 


0.6467 


0.6529 




30 


2 


0.8936 


0.8619 


0.6384 


0.6447 




30 


3 


0.8664 


0.8423 


0.6258 


0.6322 




30 


4 


0.8401 


0.8181 


0.6247 


0.6313 




30 


5 


0.8145 


0.8013 


0.6179 


0.6179 




30 


6 


0.7894 


0.7833 


0.6075 


0.6142 




20 


0 


0.9540 


0.9315 


0.7452 


0.7500 




20 


1 


0.9222 


0.9083 


0.7366 


0.7416 




20 


2 


0.8936 


0.8864 


0.7277 


0.7328 




20 


3 


0.8664 


0.8633 


0.7192 


0.7244 




20 


4 


0.8401 


0.8414 


0.7074 


0.7128 




20 


5 


0.8145 


0.8238 


0.7012 


0.7067 




20 


6 


0.7894 


0.8040 


0.6924 


0.6981 
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901 Lower Confidence Limits on RS 



Bl(r; 


5.0,0.10) 


converts 


to B2 (r ; 


39.5,1) 






Fifty 


tests for 


each component 










Classical 




Bayesian 




k 


Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


10 


0 


0.9540 


0.9600 


0.8478 


0.8508 


Same 


10 


1 


0.9222 


0.9360 


0.8413 


0.8445 


as 

Case 3 


10 


2 


0.8936 


0.9084 


0.8306 


0.8340 




10 


3 


0.8664 


0.8832 


0.8156 


0.8193 




10 


4 


0.8401 


0.8645 


0.8069 


0.8107 




10 


5 


0.8145 


0.8456 


0.7949 


0.7990 




10 


6 


0.7894 


0.8230 


0.7822 


0.7865 




5 


0 


0.9540 


0.9730 


0.9163 


0.9181 




5 


1 


0.9222 


0.9446 


0.9023 


0.9044 




5 


2 


0.8936 


0.9207 


0.8850 


0.8873 




5 


3 


0.8664 


0.9127 


0.8763 


0.8741 




5 


4 


0.8401 


0.8715 


0.8637 


0.8665 




5 


5 


0.8145 


0.8531 


0.8466 


0.8497 




5 


6 


0.7894 


0.8304 


0.8400 


0.8433 





30 



90% Lower Confidence Limits on RS 



B1 (r ; 5 . 0 , 0 . 10) converts to B2(r;39.5,l) 
Seventy-five tests for each component 

Classical Bayesian 



k 


Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


40 


0 


0.9693 


0.9193 


0.6551 


Same 


Same 












as 


as 


40 


1 


0.9481 


0.9041 


0.6541 


Case 2 


Case 2 


40 


2 


0.9290 


0.8918 


0.6419 






40 


3 


0.9109 


0.8735 


0.6369 






40 


4 


0.8934 


0.8641 


0.6326 






40 


5 


0.8763 


0.8520 


0.6271 






40 


6 


0 . 8596 


0.8365 


0.6190 






30 


0 


0.9693 


0.9372 


0.7234 






30 


1 


0.9481 


0.9201 


0.7166 






30 


2 


0.9290 


0.9044 


0.7096 






30 


3 


0.9101 


0.8884 


0.6990 






30 


4 


0.8934 


0.8782 


0.6981 






30 


5 


0.8763 


0.8626 


0.6922 






30 


6 


0.8596 


0.8509 


0.6834 






20 


0 


0.9693 


0.9541 


0.7986 






20 


1 


0.9481 


0.9382 


0.7916 






20 


2 


0.9290 


0.9221 


0.7843 






20 


3 


0.9109 


0.9033 


0.7773 






20 


4 


0.8934 


0.8925 


0.7677 






20 


5 


0.8763 


0.8800 


0.7625 






20 


6 


0.8596 


0.8616 


0.7552 
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90% Lcwer Confidence Limits on RS 



B1 (r ; 5. 0, 0. 10) converts to B2(r;39.5,l) 
Seventy-five tests for each component 

Classical Bayesian 



k 


Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


10 


0 


0.9693 


0.9718 


0.8814 


Same 


Same 












as 


as 


10 


1 


0.9481 


0.9520 


0.8763 


Case 2 


Case 2 


10 


2 


0.9290 


0.9369 


0.8677 






10 


3 


0.9109 


0.9167 


0.8559 






10 


4 


0.8934 


0.9057 


0.8488 






10 


5 


0.8763 


0.8912 


0.8393 






10 


6 


0.8596 


0.8758 


0.8290 






5 


0 


0.9693 


0.9837 


0.9354 






5 


1 


0.9481 


0.9604 


0.9244 






5 


2 


0.9290 


0.9446 


0.9108 






5 


3 


0.9190 


0.9257 


0.9042 






5 


4 


0.8934 


0.9122 


0.8942 






5 


5 


0.8763 


0.8989 


0.8806 






5 


6 


0.8596 


0.8825 


0.8756 
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90% Lower Confidence Limits on RS 



B1 (r ; 5 . 0 , 0 . 10) converts to B2(r;39.5,l) 

One hundred tests for each component 

Classical Bayesian 



k 


Failures 


Case 


Cas e 1 


Case 2 


Case 3 


Case 4 


40 


0 


0.9770 


0.9411 


0.7067 


Same 


Same 












as 


as 


40 


1 


0.9611 


0.9266 


0.7058 


Case 2 


Case 2 


40 


2 


0.9458 


0.9130 


0.6950 






40 


3 


0.9332 


0.9045 


0.6905 






40 


4 


0.9200 


0.8925 


0.6867 






40 


5 


0.9073 


0.8855 


0.6819 






40 


6 


0.8947 


0.8743 


0.6747 






30 


0 


0.9770 


0.9495 


0.7666 






30 


1 


0.9611 


0.9372 


0.7608 






30 


2 


0.9458 


0.9258 


0.7547 






30 


3 


0.9332 


0.9164 


0.7454 






30 


4 


0.9200 


0.9034 


0.7446 






30 


5 


0.9073 


0.8960 


0.7394 






30 


6 


0.8947 


0.8838 


0.7317 






20 


0 


0.9770 


0.9652 


0.8315 






20 


1 


0.9611 


0.9517 


0.8255 






20 


2 


0.9458 


0.9379 


0.8193 






20 


3 


0.9332 


0.9276 


0.8132 






20 


4 


0.9200 


0.9175 


0.8050 






20 


5 


0.9073 


0.9074 


0.8006 






20 


6 


0.8947 


0.8969 


0.7943 
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B1 (r ; 5 . 0 , 0 . 10) converts to B2(r;39.5,l) 

One hundred tests for each component 

Classical Bayesian 



k 


Failures 


Case 


Case 1 


Case 2 


Case 3 


Case 4 


10 


0 


0.9770 


0.9797 


0.9016 


Same 


Same 












as 


as 


10 


1 


0.9611 


0.9619 


0.8973 


Case 2 


Case 2 


10 


2 


0.9458 


0.9505 


0.8901 






10 


3 


0.9332 


0.9409 


0.8802 






10 


4 


0.9200 


0.9280 


0.8742 






10 


5 


0.9073 


0.9148 


0.8661 






10 


6 


0.8947 


0.9062 


0.8575 






5 


0 


0.9770 


0.9867 


0.9466 






5 


1 


0.9611 


0.9696 


0.9376 






5 


2 


0.9458 


0.9560 


0.9262 






5 


3 


0.9332 


0.9450 


0.9207 






5 


4 


0.9200 


0.9319 


0.9124 






5 


5 


0.9073 


0.9217 


0.9010 






5 


6 


0.8947 


0.9112 


0.8969 
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V. CONCLUSIONS 



The results of the simulation indicate that the classical 
approach to system reliability usually yields higher values 
of lower confidence limits than the Bayesian approach. When 
there are 10 or more components, only the beta priors with 
b parameters less than one provided higher values of lower 
confidence limits than the classical approach. When using 
the Bayesian approach, the effect of each prior density 
assumption is to obtain a more optimistic posterior density 
on each component reliability. However, there is some proba- 
bility mass that is still assigned to small values of each 
. Thus, as the number of components increase the probabil- 
ity of small values of RS increases since RS is a product of 
the R^ ’ s . Therefore, as the number of components increase, 
the effect of the prior density assumptions tends to diminish. 
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APPENDIX A 



BETA DENSITY CONVERSIONS 

The following table lists values of a*, where a* is such 
that the 10th percentile point of a B(x;a,b) distribution is 
the same as the 10th percentile point of a B(x;a*,l) distri- 
bution. 



Table of a* ' s 



b 


a=l 


a=5 


o 
1 — 1 
II 

aJ 


1.00 


0.0 


5.0 


10.0 


0.90 


1.05 


5.36 


10.77 


0.80 


1.10 


5.81 


11.72 


0.70 


1.17 


6.37 


12.91 


0.60 


1.26 


7.09 


14.46 


0.50 


1.39 


8.11 


16.61 


0.45 


1.47 


8.77 


18.02 


0.40 


1.57 


9.61 


19.78 


0.35 


1.71 


10.68 


22.06 


0.30 


1.89 


12.14 


25.13 


0.25 


2.16 


14.23 


29.55 


0.20 


2.58 


17.56 


36.57 


0.18 


2.82 


19.52 


40.70 


0.16 


3.16 


22.11 


46.16 


0.15 


3.37 


23.75 


49.59 


0.14 


3.61 


25.69 


53.69 


0.12 


4.28 


30.96 


64.78 
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Table of a*'s (Continued) 



b 


a=l 


a=5 


a=10 


0.10 


5.37 


39.46 


82.67 


0.09 


6.20 


45.98 


96.38 


0.08 


7.38 


55.24 


115.87 


0.07 


9.17 


69.29 


145.44 


0.06 


12.14 


92.61 


194.51 


0.05 


17.76 


136.82 


287.54 


0.04 


30.91 


240.42 


505.53 


0.03 


76.02 


596.85 


1256.03 


0.02 


445.64 


3528.06 


7438.38 
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COMPUTER PROGRAM 



MAIN PROGRAM 

THE PURPOSE OF THIS PROGRAM IS TO SIMULATE A SERIES 
SYSTEM OF K COMPONENTS AND TO COMPUTE BAYESIAN LOWER 
CONFIDENCE LIMITS OF SYSTEM RELIABILITY. 

INPUTS — NUM3ER OF COMPONENTS 
— BETA PRIORS 
— TEST DATA 



DOUBLE PRECISION DLBET At XX , PP » AA 

DIMENSION RSI ( 500) ,RS2 1(500) ,RS22(500) , RS23(500) , 
1NI(50),NSI(50),BETA1(50,2),3ETA2(50,2),3ETA21(50,2), 
2BETA22( 5 0, 2) , 3 E T A23 ( 5 0 , 2 ) , BET All (50,2) ,XB(40,500) , 

3 BET A3 1( 50, 2) , 3ETA32( 50, 2) , BETA 33 ( 50,2) 

DATA RSI/500 * 1.0/ 

INITIALIZE PROGRAM 

H = 26.4811 
DD1 = 47.7160 
DD = 19.8219 
IRUN = 0 
KTIME = 0 
PWANT = 0.10 

IX = 36941 
IFIRST = IX 

READ NUMBER OF COMPONENTS IN SYSTEM 

READ (5,9000) K 

9000 FORMAT ( I 3) 

READ INITIAL BETA PRIORS 

READ (5,9001) ( (BETA1( I,J) ,J=1,2) , I=1,K) 

9001 FORMAT (2F10.4) 

READ TEST DATA— NUMBER OF TRIALS AND SUCCESSES 

READ (5,9002) UNI (I ) ,NSI (I) ,1=1, K) 

9002 FORMAT (2110) 

CONVERT INITIAL BETA PRIORS TO BET A( A , 1 ) DENSITIES 

1000 CONTINUE 

DO 1099 I = 1 , K 
B = BETA1( 1,2) 

A = BETA1 (1,1) 

IF(I.EQ.l) GO TO 1011 

IF( (A.EQ.BETA1 ( 1-1, 1) ) .AND. (B.EQ.BETAK 1-1,2) ) ) GO TO 
11050 

10 11 CONTINUE 

X = 0.0 
1015 XI = 0.1 

I FLAG = 0 
X = X + X 1 

IF(X.GT.l.O) SO TO 1030 

1019 IP = 1 
IER = 1 

1020 CALL 80TR(X,A,B,P,D,IER,IP,DLBETA) 

IF (P.GT. PWANT) GO TO 1030 

x = X + XI 

IF(X.GT.l.O) SO TO 1030 
GO TO 1019 
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1030 

1040 

1050 

1099 



1100 

1101 

11 10 
1111 



11 14 



1115 

1120 



I FLAG = I F LAG + 1 

IF ( I F L AG . EQ . 6 ) GO TO 1040 

X = X - XI 

XI = XI * 0.1 
GO TO 1019 
CONTINUE 

PP = DBLE(P) 

XX = DBLE(X) 

AA = DL0G(PP)/DL0G(XX) 

A4B 1 = SNGL (AA) 

BET A2 (III) = A 4B 1 



BETA2 ( I ,2) 
GO TO 1099 
BET A2 ( 1,1) 
BETA2( I ,2) 
CONT INUE 



= 1.0 

= BETA2( 1-1,1) 
= 1.0 



ADJUST CONVERTED BETA PRIORS 

DO 1100 1=1, K 
BETA2K 1,2) =1.0 
BETA22( 1,2) = 1.0 
BETA23 (1,2) = 1.0 
BETA2K I ,1) = 0.75 * NI (I) 

BcTA22( 1,1) =1.0 * NI ( I ) 

BETA23 (1,1) = 1.5 - N I ( I ) 

IF(BETA21( I , 1) .GT. BETA2( I ,1) ) BETA21 ( I ,1 ) 
IF(BETA22( 1,1) .GT. BET A 2(1,1) ) BE TA22 ( 1,1) 
IF (BETA23( I , 1) . GT. BET A2 ( I , 1 ) ) BET A23 (1,1) 
CONTINUE 



= BETA2 (1,1) 
= 8ETA2( 1,1) 
= BETA2( I , 1) 



COMPUTE BETA PORTERIORS FROM ADJUSTED PRIORS 



CONT INUE 
DO 1110 I = 1 , K 
BETAS 1( 1,1) 

BET A32 ( I ,1 ) 

BE TA 33 ( 1,1) 

BET A3 1 ( 1,2) 
BETA32 ( I ,2 ) 

BET A33 ( I ,2) 

CONT INUE 



6ETA2K 1,1) + NSI ( I ) 

B ETA2 2 ( 1,1) +- NS I ( I ) 

8ETA23( I ,1) + NSI ( I ) 

= BETA2K 1,2) + NI ( I ) - NSI ( I ) 

= BETA2 2 ( 1,2) + NI { I ) - NS I ( I ) 

= 3 E TA2 3 ( I ,2) + NI ( I ) - NSI ( I ) 



GENERATE RANDOM BETA VARIATES AND REALIZE VALUES OF RS 

CONTINUE 
I X = IE IRST 
DO 1125 K 1 = 1,500 
RSLN1 = 0.0 
RSLN2 = 0.0 
RSLN3 = 0.0 
DO 1120 I = 1 , K 
I B 1 = BETA 3 1 ( I ,2) 

B1 = BET A3 1 ( 1,1) 

B2 = BE TA32 ( I , 1) 

B3 = B ET A33 ( 1,1) 

DO 1115 J = 1, IB1 
CALL RANDU ( IX, I Y , RAN ) 

IX = IY 

IF(RAN.EQ.O.O) GO TO 1114 

X = -A LOG (RAN) 

X.l = X/ ( J-l+Bl ) 

X2 = X/ ( J-1 + B2 ) 

X3 = X/U-1 + B3) 

RSLN1 = RSLN1 + XI 
RSLN2 = RSLN2 + X2 
RSLN3 = RSLN3 + X3 
CONTINUE 



CONT INUE 

RS21 (K 1 ) 
RS22 ( K1 ) 
RS23 ( K 1 ) 



EXP (-RSLN1 ) 
EXP (-RSLN2 ) 
EXP (-RSLN3) 
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1125 CONTINUE 

NPASS = 499 

DO 1150 I = 1, NPASS 

NSTOP = NPASS - I + 1 

DO 1150 J = 1, NSTOP 

IF(RS21(J) .LE.RS2KJ + 1 )) GO TO 1130 

TEMP = RS21(J) 

RS21(J) = RS211J+1) 

RS2KJ + 1) = TEMP 
1130 CONTINUE 

IFIRS221J) .LE.RS221 J+l) ) GO TO 1135 
TEMP = RS22 ( J ) 

RS22 ( J ) = RS22 ( J+ 1 ) 

RS22 ( J + 1 ) = TEMP 
1135 CONTINUE 

IF(RS23( J) .LE. RS23( J+l) ) GO TO 1150 
TEMP = RS23 ( J ) 

RS23 ( J ) = RS23 (J+l ) 

RS23U+1) = TEMP 
1150 CONTINUE 
2000 CONTINUE 

IF (K.LT.40) 30 TO 2201 
IR = IFIRST 

COMPUTE BETA POSTERIORS USING THE INITAL BETA PRIORS 
DO 2010 I = It K 

BETA 1 1( I t 1 ) = BETAl(Itl) + NSI(I) 

BET A 1 1 ( I 1 2 ) = BET A 1 ( 1 1 2 ) + Midi - NSI(I) 

2010 CONTINUE 

GENERATE RANDOM BETA VARIATES AND REALIZE VALUES OF RS 

I F ( I RUN . EQ . 0 ) KRUN = K 
IF CIRUN.EQ.l ) KRUN = 1 
I F ( IRUN.GT. 1 ) KRUN = 2 
I F { IRUN.EQ.5) KRUN = 3 
I F ( I RUN. EQ. 6) KRUN = 4 
DO 2200 I = It KRUN 
X = 0.999 
A = BETA 1 1 ( I ,1 ) 

B = BETA 1 1 ( 1,2) 

KK = 0 

IFtB.GE.2. 0) 30 TO 2020 
IF(B.GE.l.O) < K = 1 
GO TO 2025 
2020 KK = 2 
2025 CONT INUE 
I P = 1 
I ER = 1 

CALL BDTR ( X, A, B , P, D, I E R , 1 P , DLB E TA ) 

IP = 0 
U = P 

IF(KK.EQ.O) HT = H 
I F ( KK . EQ.l ) HT = DD1 
IFIKK.EQ.2) HT = DD 
X = 0.9 
IER = 1 

CALL BDTR(X,A,BtPtD,IER,IP,DLBETA) 

H2 = D 
X = 0.99 
IER = 1 

CALL BOTR ( X, A, B,P, D, IER, IP ,DLBETA) 

H3 = D 

DO 2120 K 1 = 1 , 500 
CALL RANDU ( IX, IY, RANI ) 

I X = I Y 

IF (RANI - U) 1, 1, 60 
1 CALL RANDU ( I X, I Y, RANI) 

IX = IY 

IF (RANI - 0.999) 6,6,1 
6 IF (RANI - 0.9) 20,20,30 
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10 CALL RANOUt IX, IY,RAN1) 

IX = IY 

IF (RANI - 0.999) 11,11,10 

11 IF (RANI - 0.9) 20,20,30 
20 CALL RANOUt IX, IY,RAN2) 

I X = I Y 

Y = RANI * (H2/0.9) 

IF ( ( RAN2 * HT)- Y) 40,10,10 

30 CALL RANDUdX, IY,RAN2) 

IX = IY 

I F( KK . EQ.2 ) GO TO 40 
I F ( RAN 1 - 0.99) 31 ,31,40 

31 CONTINUE 

Y = (RANI - 0.9) 4 ( ( H3 - H2J/0.09) + H2 

IF ( ( P.AN2 4 HT)- Y) 40,10,10 

40 X = RANI 

CALL BDTR( X,A,6,P,D,IER,IP,DLBETA) 

HI = 0 

I F ( ( RAN 2 * HT) - HI ) 5 0, 1, 1 
50 XI = RANI 

XB ( I , K 1 ) = XI 
GO TO 2120 
60 XI = 1.0 

XB ( I , K 1 ) = XI 
2120 CONTINUE 

2200 CONTINUE 

2201 CONTINUE 

DO 2250 I = 1,500 

RS 1 ( I ) = 1.0 

DO 2240 J = 1, K 

RSI ( I ) = R S 1 ( I ) 4 XB( J, I) 

2240 CONTINUE 
2250 CONTINUE 

NPASS = 499 

DO 2300 L = 1, NPASS 

NSTOP = NPASS - L + 1 

DO 2300 J = 1, NSTOP 

IF(RS1(J).LE.RS1(J+1)) GO TO 2300 

TEMP = RSI ( J) 

RS 1 ( J ) = RSKJ + l) 

RSI ( J + l ) = TEMP 
2300 CONTINUE 

WRITE RESULTS OF SIMULATION 
WR I T E ( 6 , 9300 ) 

9300 FORMAT (• 1* ,* COMPONENT BETA PRIOR CONVERTED TO B ( A , ' , 
1 • 1) TESTS SUCCESSES' ,/) 

DO 2900 I = 1 , K 

WRITE (6, 93 01 ) I , BETA! ( I , 1 ) , BET A1 ( I ,2 ) , BETA2( I , 1 ) , 
1BETA2( I ,2) .NIC I ) * N S I C I ) 

9301 FORMAT (5X, 12 B ( • , F4 . 1 , ' , ' , F 4 . 2, ' ) B ( ' , F 6 . 1 , • , ' , 

1F3. 1 , 1 ) ' , I 3 , 8X , I 3 , / ) 

2900 CONTINUE 
3000 CONTINUE 
92 00 FORMAT Cl') 

3100 CONTINUE 

9010 FORMAT (' O' ,' CASE ',11,/, 

1' 10TH PERCENTILE POINT OF RS = ',F7.4,//, 

2' 20TH PERCENTILE POINT OF RS = ',F7.4,//J 

II = 2 

WRITE ( 6,9010) I I ,RS21 ( 50) ,RS21 ( 100) 

11 = 3 

WRITE (6, 90 10) I I ,RS22(50) , RS22 ( 100 ) 

II = 4 

WRITE (6,90 10) I I ,RS23( 50) ,RS23( 100) 

RESET PARAMETERS 

KTIME = KTIME + 1 
IF(KTIME.EQ.l) K=30 
IF( KTIME. EQ.2) K=20 



41 







V 






oooooooooooooonnoooooooooooooooooooooooooooooooooooon 



IF{KTIME.EQ.3) K=10 
I F( KT IME .EQ. 4) K = 5 
I F ( KT I ME . E 3. 5 ) GO TO 9950 
GO TO 1111 
9950 CONTINUE 

IRUN = IRUN + 1 
IF( IRUN.EQ.7) GO TO 9999 
RE AD (5, 9002 ) ( M I ( I) , NS I ( I ) , I = 1, 5 ) 
K — 40 
KT I ME = 0 
GO TO 1101 
9999 CONTINUE 

9998 FORMAT ('1’,'RUN COMPLETE',/) 
WRITE (6,9998) 

STOP 
END 

SUBROUTINE BDTR 



PURPOSE 

COMPUTES FIX) = PROBABILITY THAT THE RANDOM VARIBLE 
U, DISTRIBJTED ACCORDING TO THE BETA DISTRIBUTION 
WITH PARAMETERS A AND B, IS LESS THAN OR EQUAL TO 
X. F ( A , B , X ) , THE CRCINATE OF THE BETA DENSITY AT 
X IS ALSO COMPUTED. 

USAGE 

CALL BDTR(X,A,B,P,D, IER) 

DISCRIPTION OF PARAMETERS 

X - INPUT SCALAR FOR WHICH P ( X) IS COMPUTED. 

A - BETA DISTRIBUTION PARAMETER (CONTINUOUS). 

B - BETA DISTRIBUTION PARAMETER (CONTINUOUS). 

P - OUT P J T PROBABILITY. 

D - OUTPUT DENSITY. 

IER - RESULTANT ERROR CODE WHERE 

IER= 0 NO ERROR 

I ER=-1 , +1 CDTR HAS BEEN CALLED AND AN ERROR 
HAS OCCURRED. SEE CDTR. 

I ER=-2 AN INPUT PARAMETER IS INVALID. 

I ER=+2 INVALID OUTPUT. 

REMARKS 

SEE MATEMATICAL DESCRIPTION. 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
DL GAM 
NOT R 
CDTR 

METHOD 

REFER TO R. E. BARGMANN AND S.P. GHOSH, STATISTICAL 
DISTRIBUTION PROGRAMS FOR A COMPUTER LANGUAGE, 

IBM RESEARCH REPOST RC-1094, 1963. 



MODIFICATION TO SUBROUTINE BDTR 

USAGE CALL BDTR ( X , A , B , P , C , I ER , I P , D LBE TA ) 

DISCRIPTION OF ADDITIONAL PARAMETERS 

IP ACTS AS A FLAG TO ELIMINATE THE 

COMPUTATION CF DLBETA. 

DLBETA — ADDED TO PARAMETER LIST TO SAVE. 

IER UTILIZED TO EXIT SUBROUTINE AFTER 

DENSITY IS COMPUTED IF P(X) IS NOT 
REQUI RED. 
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SUBROUTINE BDTR ( X, A , B , P , D, IE R, I P , DLB ET A ) 

DOUBLE PRECISION XX , DL XX , D L 1 X , A A , BB , G1 , G2 , G3 , G4 , DD, P P , 
1X0,FF,FN,XI,SS,CC,RR,DLBETA 

TEST FOR VALID INPUT DATA 

MODIFICATION — THE FOLLOWING CARD WAS OMITcD. 

I F C A-( .5-1. E-5) ) 640,10,10 

10 CONTINUE 

20 I F { A- 1 . E + 5 ) 30,30,640 
30 IFtB-l.E+5) 40,40,640 
40 I F ( X ) 640,50,50 
50 IF(l.-X) 640,60,60 

COMPUTE LOG { B ET A { A , B ) ) 

60 AA=DBL E ( A ) 

BB=DBLE ( B ) 

MODIFICATION — THE FOLLOWING CARD WAS ADDED. 

IF ( IP .EQ.O) GO TO 65 

CALL DLGAM ( AA , G1 ,I0K) 

CALL DLGAMIBB, G2, IOK) 

CALL DLGAMI AA+3B,G3, ICK) 

DLBETA=G1+G2-G3 

MODIFICATION — THE FOLLOWING CARD WAS ADDED. 

65 CONTINUE 

TEST FOR X NEAR 0.0 OR 1.0 

IFIX-l.E-8) 80,80,70 
70 IFl ( l.-XJ — 1.E-8J 130,130,140 
80 P=0.0 

IFIA-l.) 90,100,120 
90 D= 1 . E + 7 5 
GO TO 660 
100 DD=-DLBETA 

IFIDD+l .68D02) 120,120,110 

110 DD=DEXP { DD ) 

D=SNGL ( DO ) 

GO TO 660 
120 D=0 . 0 

GO TO 660 
130 P=1.0 

1FIB-1.) 90,100,120 

SET PROGRAM PARAMETERS 

140 XX=DBL E ( X) 

DL XX=DL0G ( XX ) 

DL1X = DL0G( 1 .DO-XX ) 

XO=XX/( 1. DO-XX) 

ID=0 

COMPUTE OROINATE 

DD = ( A A-l . DO ) *D LXX+ ( B 3- 1 . DO ) DL IX- DL B ETA 
I FIDO- 1. 68D02) 150,150,160 

150 I F ( DD + 1 . 68 DO 2 ) 170,170, 180 

160 D=1 . E75 

GO TO 190 
170 D=0 . 0 

GO TO 190 
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180 DD=D£XP(DD) 

D=SNGL COD) 

MODIFICATION — THE FOLLOWING CARD WAS ADDED. 

IF ( IER. EQ.O) GO TO 670 

A OR B OR BOTH WITHIN l.E -8 OF 1.0 

190 I F ( AB S ( A-l. )-l . E— 8 ) 200,200,210 

200 IF (ABS( a-l. ) -1 . E- 8 ) 220,220,230 

210 I F { ABS( B-l.J-l.E- 8 ) 260,260,290 

220 P=X 

GO TO 660 
230 PP=B3*DL1X 

I F ( PP + 1 .08 D02 ) 240,240,250 
240 P=l. 0 

GO TO 660 
250 PP=DEX P ( PP ) 

PP = 1 . D 0-PP 
P=SNGL ( PP) 

GO TO 600 
260 PP=AA*DLXX 

I F { PP+1 . 68 D02 ) 270,270, 280 

270 P=0.0 

GO TO 660 
280 PP = DEX P ( PP ) 

P=SNGL ( PP) 

GO TO 600 

TEST FOR A OR B GREATER THAN 1000.0 

290 I F ( A- 1000 . ) 300,300,310 
300 IF ( B-iOOO. ) 330,330,320 

310 XX=2. DO*AA /XO 
XS=SNGL ( XX ) 

A A=2. DO*BB 
DF=SNGL ( AA ) 

CALL CDTR( XS, DF, P, DUM'Y , IER) 

P=l. 0-P 
GO TO 670 
320 X X = 2 . DO ^ B B ‘~X 0 
XS=SNGL( XX) 

A A=2 . DO X ’ : AA 
DF=SNGL ( AA ) 

CALL CD TR( XS , 0 F , P , DUMMY , I ER) 

GO TO 670 

SELECT PARAMETERS FOR CONTINUED FRACTION COMPUTATION 

3 30 I F ( X— . 5 > 340,340,380 
340 IF(AA-l.DO) 350,350,360 
350 RR = AA+1 . DO 
GO TO 370 
360 RR=AA 
370 DD=DLXX/5.D0 
DD=DEXP ( DD) 

DD= ( RR- l..DO)-( RR+BB-1. DO)-XX*DD +2. DO 
IF(DD) 420,420,430 
380 IF(BB-l.DO) 390,390,400 
390 RR=B3+ 1 . DO 
GO TO 410 
400 RR = B B 
410 DD= DL 1 X / 5 .DO 
DD=DEXP { DD ) 

DD = ( RR- 1 . D 0) - ( AA+RR-1 . DO) * (1 .DO-XX)*DD +2 . DO 
IF(DD) 430,430,420 
420 I D = 1 

FF=DL1X 
DL 1X=DL)'.X 
DLXX=FF 
X0=1. DO/XO 
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non non non 



FF=AA 
AA = BB 
BB=FF 
G2 = G 1 

TEST FOR A LESS THAN 1.0 
430 FF=0 .DO 

IF(AA-l.DO) 440,440,470 
440 CALL DLGAM(AA+1.D0,G4, I OK) 

DD=AA*DLXX+BB*DL1X+G3-G2-G4 
I F ( DD+- 1 . 68D02) 460,460,450 

450 FF=FF+DEXP ( DD) 

460 AA=AA+ 1 . DO 

COMPUTE P JSING CONTINUED FRACTION EXPANSION 

470 FN=AA+BB- 1 . DO 
RR=AA-1 .00 
I I =80 

XI =DFL0AT ( I I ) 

SS= ( ( BB-XI )*( RR + XI ) )/{ (RR+2.D0*XI-1.D0) * ( RR + 2 . DO* XI ) ) 

SS=SS*X0 

DO 480 1=1,79 

I 1 = 80-1 

XI =DFLCAT ( I I ) 

DD= { X I* ( FN+XI ) ) /( ( RR+2 . DO- X I + 1 . DO) * ( RR+2 . DO* XI ) ) 
DD=DD*X0 

CC=( 1 8B-XI )*( RR+XI ) )/( (RR*2. D0*XI-1. DO)* iRR+2.D0*XI ) ) 
CC=CC*X0 

SS=CC/ ( 1 .DO+DD/ { 1. DO-SS ) ) 

430 CONTINUE 

SS=1.D0/(1. DO-SS) 

IF (SS) 650,65 0,490 
490 CALL DLGAM(AA*BB,Gl,IfiK) 

CALL DLGAM ( AA+ 1 . DO , G4- I OK) 

CC=G1-G2-G4+AA*DLXX+(3B-1 .DO) *DLiX 
PP=CC+DLOG ( SS ' 

IF{PP + 1 .68 DO 2. 500,500,510 

500 PP=FF 

GO TO 520 

510 PP=DEXP ( PPJ+FF 
520 IF(ID) 540,540,530 
530 P D = 1 . DO— PP 
540 P=SNG L ( PP) 

SET ERROR INDICATOR 

I F ( P ) 550,570,570 
550 IF(A3S( P J-i.E-7) 560,560,650 
560 P=0.0 

GO TO 660 

570 IF(l.-P) 580,600,600 

530 IF (ABS(l.-P)-i .E-7) 590,590,650 

590 P=1.0 

GO TO 660 

600 IF(P-l.E-S) 610,610,620 
610 P=0.0 

GO TO 6o0 

620 I F ( (1.0-P)— l.E-8) 630,630,660 
630 P=1.0 

GO TO 660 
640 I E R=- 2 

D=-l . E75 
P=-l .E75 
GO TO 670 
650 I ER = +2 

P= 1.E75 
GO TO 670 
660 I ER = 0 
670 RETURN 
END 
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